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The first normal stress difference (A/i) and tlie microstructure in a dense sheared 
granular fluid of smooth inelastic hard-disks are probed using event-driven simula- 
tions. While the anisotropy in the second moment of fluctuation velocity, which is a 
Burnett-order effect, is known to be the progenitor of normal stress differences in di- 
lute granular fluids, we show here that the collisional anisotropics are responsible for 
the normal stress behaviour in the dense limit. As in the elastic hard-sphere fluids. 
Ml remains positive (if the stress is defined in the compressive sense) for dilute and 
moderately dense flows, but becomes negative above a critical density, depending on 
the restitution coefficient. This sign-reversal of A/i occurs due to the microstruc- 
tural reorganization of the particles, which can be correlated with a preferred value 
of the average collision angle Oav = '7r/4±7r/2 in the direction opposing the shear. 
We also report on the shear-induced crj/staZ-formation, signalling the onset of fluid- 
solid coexistence in dense granular fluids. Different approaches to take into account 
the normal stress differences are discussed in the framework of the relaxation- type 
rheological models. 
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I. INTRODUCTION 



In the last decade, a lot of research activity took place to unveil the properties of granular 
materials^'^, primarily because of their industrial importance, but also due to their fascinat- 
ing properties. This has unraveled many interesting and so far unresolved phenomena (for 
example, clustering, size-segregation, avalanches, the coexistence of gas, liquid and solid, 
etc.). Under highly excited conditions, granular materials behave as a fluid, with prominent 
non-Newtonian properties, like the normal stress differences^. While the normal stress dif- 
ferences are of infinitesimal magnitudes in a simple fluid (e.g. air and water), they can be of 
the order of its isotropic pressure in a dilute granular gas^. From the modelling viewpoint, 
the presence of such large normal-stress differences readily calls for higher-order constitutive 
models^'^ even at the minimal level. 

Studying the non-Newtonian behaviour is itself an important issue, since the normal 
stresses are known to be the progenitors of many interesting and unique flow-features (e.g. 
rod-climbing or Weissenberg-effect, die-swelling, secondary flows, etc.'^) in non-Newtonian 
fluids. Also, normal stresses can support additional instability modes (for example, in 
polymeric fluids and suspensions^" which might, in turn, explain some flow- features of 
granular fluids. For example, particle-clustering^^""^^ has recently been explained from the 
instability-viewpoint using the standard Newtonian model for the stress tensor^^'^^'^^ 

The kinetic theory of Jenkins & Richman^^ first showed that the anisotropy in the second 
moment of the fluctuation velocities, due to the inelasticity of particle collisions, is responsi- 
ble for such normal stress behaviour. They predicted that the flrst normal stress difference 
(deflned as A/i = {Uxx — ^yy)/P: where U^x and Uyy are the streamwise and the transverse 
components of the stress dcviator, respectively, and p is the isotropic pressure, see section 
JIB) is maximum in the dilute limit, decreases in magnitude with density, and eventually 
approaches zero in the dense limit. Goldhirsch & Sela^ later showed that the normal stress 
differences appear only at the Burnett-order-description of the Chapman- Enskog expansion 
of the Boltzmann equation. Their work has clearly established that the origin of this effect 
(in the dilute limit) is universal in both atomic and granular fluids, with inelasticity playing 
the role of a magnifier and thus making it a sizeable effect in granular fluids. While the 
source of the normal stress differences in the dilute limit has been elucidated both theo- 
retically and by simulation, its dense counterpart has not received similar attention so far. 
This is an important limit since the onset of dilatancy (volume expansion due to shcar^'''^^), 
crystallization, etc. occur in the dense regime, which in turn would influence the normal 
stress differences. 

Previous hard-sphere simulations^^'^'^^ did look at the normal stress differences, but they 
did not probe the dense limit in a systematic way. These simulations showed that the flrst 
normal stress difference vanishes in the dense limit, in line with the theoretical predictions of 
Jenkins & Richman^^. On the other hand, the soft-sphere simulations of Walton & Braun^*^, 
with frictional particles, showed that this quantity can change sign in the same limit. Our 
work with smooth inelastic hard-disks unequivocally demonstrates that J\fi, indeed, changes 
its sign at some critical density in the dense regime, due to the sign-change of its collisional 
component at a critical density, which depends on the value of the coefficient of restitution 
(e). More importantly, we show that the origin of A/i in the dense hmit is distinctly different 
from that in a dilute granular gas. At the microstructural-lcvcl, certain topological changes 
in the anisotropic structure of the collision- angle distribution with density are responsible 
for the observed sign- reversal of A/i. 
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We use the familiar smooth hard-disk model for an event-driven simulation^^ of the 
uniform shear flow configuration, focussing mainly on the normal stress behaviour and the 
microstructure formation as functions of the density and inelasticity. The details of the 
simulation technique and the relevant macroscopic quantities are described in section ITTl 
The simulation results on the first normal stress difference, the radial distribution function, 
the collision angle distribution and the crystalline-structure are presented in section lllll 
Possible modelling approaches to incorporate the normal stress differences are discussed in 
section |TVl In section |V| we summarize our findings, with suggestions for possible future 
work. 

II. SIMULATION METHOD 

We consider a collection of smooth inelastic hard-disks in a square box of size L under 
uniform shear flow — let x and y be the streamwise and transverse directions, respectively, 
with the origin of the coordinate-frame being positioned at the centre of the box. The 
snapshot of a typical simulation, with non-dimensional coordinates, is shown in Fig. 1(a). 
Note that the dimensional quantities are denoted by tildes, and the reference length, time 
and velocity scales for non-dimensionalization will be specified later in this section. 

Let the diameter and the mass of the particle be a and m, respectively. The pre- and 
post-colhsional particle velocities of particle 1 are denoted by Ci and c'^, respectively. Hence, 
the velocity of particle 2 relative to 1 is C21 = 62 — Ci. Let = k be the unit vector directed 
from the center of particle 2 to that of particle 1 at contact. The pre- and post-collisional 
velocities are related by the expression: 

k-c'ai = -e(k-C2i), (1) 

where e is the coefficient of normal restitution, with < e < 1; note that we restrict ourselves 
to perfectly smooth particles. The expression for the coUisional impulse is 

~ Tn 

I = m(c;-ci) = -(l + e)(k-C2i)k, (2) 

directed along k. 

A. Model system and algorithm 

The system is periodic in x-direction, i.e. a particle crossing the left/right boundary 
re-enters the system through the opposite boundary at the same vertical position with 
unchanged velocities. To impose a uniform shear rate (7 = U / L) in the y-direction,_the top 
and bottom image boxes, bounding the central box, are set in motion with velocities U /2 and 
— f/ /2, respectively, in the streamwise direction. This is the standard approach to attain the 
state of uniform shear flow (USF) by imparting momemtum transfer by shearing, originally 
introduced by Lees & Edwards^^. Overall, this system represents an extended doubly-periodic 
system where the periodicity in the transverse direction is in the local Lagrangian frame. 
In a typical simulation, the disks are initially placed randomly in the computational box, 
and the initial velocity field is composed of the uniform shear and a small Gaussian random 
part. An event-driven algorithm is then used to update the system in time, the details of 
which may be found in Alam & Luding^^'^"^. 
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To ascertain whether the system has reached the statistical steady- state, the time evohi- 
tion of the average fluctuation kinetic energy ('granular' energy, deflned in the next section) 
is monitored, see Fig. 1(6). Due to the balance between the shear work and the coUisional 
dissipation under homogeneous shear deformation, the granular energy attains a constant 
value at the steady state. Depending on the value of the coefficient of restitution and the 
number of particles, it takes about thousand collisions per particle to reach such a statisti- 
cal steady-state - the lower the value of e, the more quickly the system reaches the steady 
state and vice versa. The simulation is then allowed to run for another 15000 collisions per 
particle to gather data to calculate the macroscopic quantities. A few longer runs (30000 
collisions per particle) were also checked, with no significant change on the measured quan- 
tities. Another quantity which was simultaneously monitored, along with granular energy, 
is the linearity of the streamwise velocity profile across the Couette gap, and we found that 
the calculated shear rate (i.e. the slope of the velocity profile) fluctuated around the imposed 
shear rate by at most 1% at densities where crystallization is not evidenced. 



B. Macroscopic quantities 

With L, 7~^, jL, and m, used as the reference length, time, velocity, and mass, respec- 
tively, the relevant dimensionless quantities are: 



J, (c,u,C) = :^(c,u,C), P^-,^, T = — 
L pern cr 7 



where u is the 'hydrodynamic' velocity, C = c — u the fluctuation (peculiar) velocity of 
particles, p the material density of particles, P the stress tensor, and T the granular energy. 

The macroscopic stress, as measured in discrete particle simulations, is a byproduct of 
the particle-level mechanisms of momemtum transfer. As in the hard-core model of dense 
gases, the stress is the sum of its kinetic and coUisional components. The former arises 
from the transport of momentum as the particles move through the system carrying their 
momentum, while the latter is due to the direct interparticle collisions. The homogeneity 
of the uniform shear flow allows us to calculate the stress by averaging it over the whole 
computational box"''^'^'^^. 

The stress, deflned in the compressive sense, may be decomposed in the standard way: 

p = p'^ + P'^ = pi + n, (4) 

where p is the pressure, 11 the pressure deviator and 1 the unit tensor. From the off-diagonal 
components of the pressure deviator, we can calculate the shear viscosity which relates the 
rate of strain to the shear stress: 

p = -Iixy/'^- (5) 

For the steady uniform shear flow, thus, the dimensionless shear viscosity can also be in- 
terpreted as the shear stress due to our adopted scaling, du/dy = 7 = 1. The diagonal 
components of the pressure deviator can be non-zero, giving rise to normal stress differ- 
ences. The flrst normal stress difference is deflned as 

M = (^---^yy\ (6) 
p 
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Note that we have scaled this quantity by pressure to ascertain its relative magnitude with 
respect to pressure. For a standard Newtonian fluid, A/i = and thus A/i is an indicator of 
the non-Newtonian character of the fluid. Ni can be decomposed into kinetic and coUisional 
parts: 

j\f^=M^+M^=^-^ — yyl^^-^ — yyl. (7) 

p p 

Note that the sign of Mi crucially depends on the convention used to define the stress 
tensor. For example, in the rheology literature, stress is typically defined in the tensile 
sense^. A positive A/i for the compressive case is equivalent to its negative value for the 
tensile case and vice versa. This point should be kept in mind while making any comparison 
with data in the rheology literature. 

From the trace of the kinetic stress tensor, P'^, one can calculate the granular energy. 



T 



2(72 



1 ^ 



N 
i=i 



(8) 



which is a measure of the random motion of the particles with respect to the mean motion. 

There are two dimensionless control parameters: the volume fraction of particles (z/) and 
the coefficient of normal restitution (e). The shear rate is also a control parameter, however, 
due to normalization we have 7 = 7/7=1, and changing the value of 7 does not influence 
the reported results; the imposed shear rate 7 is thus kept fixed at unity. The simulations 
are carried out for the whole range of solid volume fractions, varying from the dilute to the 
dense limit, over a large range of values for the coefficient of restitution (e = 0.3-0.99). For 
most of the simulations, the number of particles are fixed to = 1024, and increasing the 
value of N by fourfold {N — 4096) did not affect the reported quantities noticeably; for 
example, the change in A/i was about 4.3% and 5.1% at z/ = 0.3 and 0.75, respectively, for 
a restitution coefficient of e = 0.9. We note here that the system-size dependence of the 
rheological quantities (pressure and viscosity) is known to be strong only for a small number 
of particles {N < lOO)^^'^^'^^ 

For the typical simulation in Fig. 1(a), at steady-state, after 2 x 10^ collisions, the 
parameter values were z/ = 0.5, = 1024 and e = 0.7. The variations of the granular 
energy T and the calculated shear rate 7ca/ with time are shown in Fig. 1(6), along with 
corresponding initial variations in two insets. (The data represent the instantaneous values 
of T and ^cai sampled at a regular interval of 400 collisions - no time averaging is involved 
here.) Note that 7^0/ was computed by binning the system into 20 equal-size bins in the 
transverse direction and taking averages over all particles in each bin. It is observed that 
the granular energy reaches its steady value (T = 0.6121 ± 0.022) quickly after the initial 
transients and the calculated shear rate fluctuates around its imposed value (7 = 1) by 
about ±1%. The fluctuations in both T and ^cai at steady-state are due to the finite-size of 
the system and diminish with increasing number of particles as A""^/^. 



III. RESULTS 



For detailed results on the quantities pressure, shear viscosity and granular energy, and 
for their comparison with kinetic theory predictions, we refer to our recent study Here, 
we will mainly focus on the behaviour of the first normal stress difference and its kinetic 
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and coUisional components. Wc also present results on the pair distribution function and 
the collision angle distribution to characterize microstructures. Lastly we will present re- 
sults on crystaZ-formation at high densities, signalling the coexistence of fluid and solid, 
complementing recent results in non-sheared systems^^'^^. 

A. Normal stress difference 

Figure 2(a) shows the variation of the first normal stress difference (A/i) with density 
for two values of the coefficient of restitution. It is observed that Mi is maximum at the 
dilute limit and decreases thereafter with u. The overall variation of Mi with density looks 
similar at other values of the coefficient of restitution, with a difference in the magnitude 
oi Ml- The inset in Fig. 2(a) shows that Mi decreases quite sharply in the dense limit and 
becomes negative at some density {v = V). Increasing the value of the restitution coefficient 
decreases this critical density V. The arrows on the left-ordinate indicate the asymptotic 
values of Mi for a two-dimensional granular gas in the dilute limit^: 

TVi = 1.0448(1 -e^). (9) 

The anisotropy in the second moment of the fluctuation velocity is primarily responsible 
for the finite normal stress difference in the dilute limit^^'^^ and this shows up only at the 
Burnett-order of the Chapman- Enskog expansion^. We should mention here that the limit 
e ^ 1 is singular and the normal stress difference survives even in the elastic limit as pointed 
out by Goldhirsch & Sela^. The corresponding expression for Mi in a molecular gas is: 

Ml ^ 1.358f^, 

where I is the mean free path and {u^) is the rms of the velocity fluctuations. However, be- 
cause of its extremely small magnitude under normal conditions, the normal stress difference 
is not measurable in a molecular fluid. 

Previous hard-sphere simulations of Campbell and coworkers^'^^ are in variance with our 
result in that they found Mi as z/ — Vmax- However, the soft-sphere simulations 
of Walton & Braun^°, with frictional particles, support our observation that Mi indeed 
undergoes a sign-reversal. To better understand what is responsible for the sign-reversal of 
jVi, we look at the kinetic and coUisional components of the first normal stress difference. 
Figure 2(6) shows the variations of M^ and Mi with density at e = 0.7. We observe that 
Ml is maximum at the dilute limit and decreases monotonically to zero as p approaches the 
packing limit. Except for the dense limit, the overall behaviour of Mi represents that of 
the total normal stress difference. The coUisional component, jVf, shows a non-monotonic 
variation with density: Mi is zero in the dilute limit, increases with increasing u, remains 
almost constant for intermediate densities, and then decays sharply in the dense limit. 
Interestingly, Mi becomes negative at some critical density {u = V) beyond which the 
behaviour of Mi mirrors that of Mi (see inset) . Thus the normal stress behaviour in the 
dense regime is clearly due to the anisotropy in the coUisional stress. 

Recall that the kinetic theory of Jenkins & Richman^^'^® predicts that Mi — in the 
dense limit. The predictions of the revised Enskog theory of Santos et al.^^ are in line with 
that of Jenkins & Richman, even though their kinetic model is claimed to be valid even in 
the crystalline-phase. Since the source of normal stress differences in all these theories is 
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linked to the anisotropy in the second moment of velocity fluctuations (which vanishes as 
— T-'max)i they are unable to predict the correct behaviour of normal stresses in the dense 
limit. We would like to stress here that, as mentioned in the Introduction, many fascinating 
non-Newtonian effects^ are primarily determined by the first normal stress difference and 
its sign. For example, the rod-climbing effect can occur in a granular fluid if Mi < 0, i.e. 
only in the dense limit (assuming that the second normal stress difference is negligible). 
Furthermore, having a constitutive model which reproduces the correct sign of A/i is also 
important, since it is well known in the rheology literature that the odd sign for J\fi leads to 
the instability of the rest state^°. 

The sign-reversal of AAf can be succinctly presented as a phase-diagram in the {u, e)-plane 
by plotting the zeros of A/'f as a function of the coefficient of restitution, see Fig. 3. Below 
the solid line, A/'f is positive, and negative above it. Also plotted in this figure is the line 
for the zeros of A/i which, as expected, hes shghtly above. Thus, the first normal stress 
difference is zero along the solid line which may be called the symmetry-line. It is observed 
that decreasing the coefficient of restitution increases the critical density (V) at which AAf 
changes sign. As we approach the elastic limit, u depends strongly on the value of e. We 
further note that as e — > 1, 1/ ~ 0.62 which is well below the freezing-point density of a 2D 
hard-disk fiuid, i// ^ O.TO^^'^^. 

We need to mention here that the sign- change of A/i is not uncommon in other non- 
Newtonian fluids. For example, in non-Brownian viscous suspensions, A/i changes sign at 
high Peclet number^^'^^. However, the reason for this effect is quite different in granular 
fluids as we show below. 

B. Microstructural features 

To understand the microstructural mechanism for the origin of the first normal stress 
difference and its sign-reversal, here we probe several microstructural features of a dense 
granular fiuid. 

1. Radial distribution function 

Typical snapshots of the system in the dense regime are shown in Fig. 4 at four different 
densities with e = 0.7. Note that the density for the subplot 4(c) is u = 0.725 for which 
AAf ~ 0. Looking at the corresponding distribution of granular energies (not shown here 
for brevity), we could find signatures of clusters (group of particles) with lower energies 
surrounded by particles with higher energy. To understand the fiow-microstructures and 
their energetics at such high densities, we need to probe the pair distribution function and 
similar measures for the structure of the packing. 

Figure 5{a-d) shows the radial distribution function g{r) at four densities, with parame- 
ter values as in Fig. 4. The thin, dotted lines are data from a non-sheared, homogeneous, 
elastic systcm^^'^^, whereas the thick, solid lines represent a sheared situation with rather 
strong dissipation e = 0.7. (Note that the elastic distribution function has been measured 
under non-sheared periodic boundary conditions; also for e somewhat smaller than unity, 
the same results were obtained as long as the system remains homogeneous.) We observe 
that the weak difference at low density u — 0.6 grows with increasing density, concerning 
two aspects: 
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(i) The peak value of contact in the sheared systems is always larger than that in a homoge- 
neous system of the same density, and the difference increases strongly with density, another 
indicator for clusterings^. 

(ii) The peaks and valleys, which allow to distinguish between different lattice structures, are 
different in the sheared case when compared to the non-sheared situation. In the former case, 
peaks at rja = 1,2,3,... are observed, signalling shell-formation about any test-particle. 
In the latter case, the peaks at r/a = 1, \/3, 2, . . . indicate a crystallization transition and 
the development of a triangular lattice. 

The peaks in the sheared situation, e.g. at r/a — 2,3, become sharper as the density 
is increased, but the one at r/o" = 4 is not well-defined even at u = 0.75; higher order 
peaks are almost invisible, indicating long-range disorder due to the shearing in contrast 
to the long-range order that evolves in the non-sheared system. The comparison between 
the sheared and non-sheared cases suggest that the structure-formation is much slower in 
a sheared fluid which, in turn, implies that, as expected, the freezing-point density of the 
former would be larger than that of the latter. The splitting of the second-peak in the non- 
sheared case corresponds to the onset of freezing transition^^; with increasing density this 
splitting becomes much more prominent, with similar structural-features appearing at the 
successively higher-order peaks. For the sheared system, however, we do not observe similar 
splitting, rather we see a sharp second-peak. This also occurs for a highly inelastic sheared 
fluid at a much lower density^'^; hence, a sharp second peak is a signature of short-range 
ordering due to the dissipative particle-clustering. The higher-order peaks in the sheared 
case become prominent only if we go beyond u — 7r/4, which corresponds to the limit of 
perfect square-packing, but it is difficult to maintain homogeneous-shearing at such high 
densities, because the system splits into two parts, a dense, cold, crystalline area and a 
dilute, hot, fluid area - see below for details. 

Figure 6 shows the variations of the pressure and viscosity functions (/p = p/pT and 
ffj. — l^/po'Vr, respectively) with density at a restitution coefficient e = 0.7. We observe 
that both increase monotonically with density, much beyond u — 0.7001, and we did not 
find the hysteretic van-der- Walls loop in our pressure data upto u = 7r/4, another indication 
that the crystallization is hindered/delayed by shear. These observations, together with our 
result that in the elastic limit J\fi changes sign at a much lower density {u ~ 0.62, see Fig. 3) 
than the corresponding freezing-density, suggest that the sign-reversal of jVi is not related 
to the freezing-transition. 

2. Collision angle distribution 

Associated with the sign-reversal of the first normal stress difference is a change in the 
relative magnitudes of the normal stress components {Pxx and Pyy)- Subtle changes in the 
direction and magnitudes of the colHsional-modc of momcmtTim transfer could infiuence 
the individTial components of the stress tensor. In order to test this hypothesis, we focus 
on the collision angle distribution function, C{9), which is defined such that C{0)d6 is the 
probability of collisions occuring at an angle lying between d and 9 + d9, with the angle 9 
being measured in the anticlockwise direction from the positive x-axis (see Fig. 13). For a 
fiuid in equilibrium, all collisions are equally likely, and hence C{9) is a uniform function of 
9, i.e. C{9) = l/n ^ 0.318309. For a non-equilibrium system (e.g. shear fiow), however, 
preferred collisions are dictated by the nature of the external field, leading to an anisotropic 
distribution for C{9y^~^^. Following Savage & Jeffrey^^ and Campbell & Brennen^^, an 
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explicit expression for C{9\ u, e) can be derived for the case of uniform shear flow as detailed 
in the Appendix. Note that the angular dependence of C{6; i/, e) will be modified by both 
the density and the restitution coefficient. 

Figure 7 shows the comparison of our simulation data on C{6) with the theoretical pre- 
dictions of equation ()A9|) for two different values of the restitution coefficient at a density 
z/ = 0.6. It is observed that the probability of collisions is higher on the up stream-faces of 
the colliding particles, i.e. for 9 G [7r/2,7r] and 9 G [— 7r/2,0] (i.e. the hatched-areas on the 
test-particle in Fig. 13). This is a consequence of the imposed shear-field which compresses 
the flow-structure along the 37r/4-direction and stretches it along the 7r/4-direction. Re- 
garding the comparison with theory, there is, clearly, a phase-difference between theory and 
simulation, and the overall agreement is only qualitative. 

Note in Fig. 7 that the probability of collisions on the upstream-faces increases further 
as the restitution coefficient increases. This, in turn, suggests that particle motion be- 
comes more streamlined (i.e. along the streamwise direction) with increasing dissipation 
levels, which will naturally lead to a reduction in the transverse component of the fluctu- 
ation velocities of the particles. Thus, the macroscopic manifestation of such microscopic 
streamlined-motion would be an increase in the magnitude of the kinetic component of the 
flrst normal stress difference (A/'f). 

Turning our attention to the range of densities where TVi undergoes a sign-reversal, we 
show the collision angle distributions C{9) in Fig. 8 as polar plots with e = 0.7; the 
corresponding densities are as in the subplots of Fig. 4. It is observed that the anisotropic 
structure of C{9) gets further modifled in this regime, with distinct peaks appearing near 9 = 
and 27r/3 (see subplot b). While the peak at = corresponds to head-on collisions between 
particles in the same-layer, the one at ^ = 27r/3 clearly signals the onset of triangular- 
structure formation. Another noteworthy point is that the collisions on the downstream-faces 
of the colliding particles are rare at these densities and hence C [9) can be approximated solely 
by its contributions from the second- and fourth- quadrants {9 G [7r/2,7r] and 9 G [— 7r/2,0], 
respectively). 

Since the momentum transfer occurs mainly due to collisions in the dense regime, the 
stress tensor can be approximated by 




(k®k)C(^)d^, 



where k is the unit vector joining the line of centres of the two colliding disks. Assuming 
now that all the collisions would occur at some average collision angle 9av so that C{9) = 
C{9av)i and recalling that C{9) is well represented in this regime by restricting 9 only in the 
second- and fourth-quadrants, the expression for the first normal stress difference simplifies 
significantly to 

ATi ~ [(k ® k). - (k ® k),],^,^^ C{9a,). (10) 

It is trivial to check that A/i = at 9av = — vr/4. From our simulation data, we have 
calculated 9av by averaging C{9) over the second- and fourth-quadrants, whose variation 
with density is plotted in Fig. 9 for two restitution coefficients. It is observed that 9av crosses 
through — 7r/4 (i.e. 37r/4) at around the critical density u for all restitution coefficients. For 
example, 9av ~ —45.16° and —45.04° at 17 = 0.725 and 0.67, respectively, where A/i changes 
sign. Thus, the microstructural signature of the sign-reversal of A/i is directly correlated 
with the average collision angle being greater or less than — 7r/4 (or Sn/A). 
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C. Crystallization: Fluid-solid coexistence 



Figure 10(a) shows a snapshot of the system at the steady state with parameter values 
being set to u — 0.8 and e — 0.9. It is observed that a sohd-layer coexists with two fluidized 
zones on either side of it. A closer look into the solid- layer reveals that the particles are 
arranged in a triangular-packing, representing a crystal, and thus we have a clear evidence 
for fluid-solid coexistence. The corresponding instantaneous streamwise velocity profile at 
t — 390 (i.e. the image-boxes have moved 390 strain units from their original position) is 
shown in Fig. 10(6); the coarse-graining is done by binning the system into 20 equal-size 
bins in the transverse direction and then taking averages over all the particles in each bin. 
Clearly, the shearing is inhomogeneous across the Couette-gap: the crystal is aligned along 
the streamwise direction and hence we call it a layered-cijstal; the shear-rate in the fluidized 
regimes on either side of the crystal is almost uniform. Note that the asymmetric nature of 
the velocity profile also signals the breakdown of the Lees-Edwards boundary condition as a 
motor for the homogeneous shear. The formation and the time-evolution of this crystal can 
be ascertained from Fig. 10(c) which shows the corresponding evolution of the streamwise 
velocity at early times. We observe that the crystal has fully formed ai t — 76, and the 
velocity profile remains antisymmetric about y — till t — 150. With further time-evolution, 
however, the crystal does not remain stationary in the transverse direction, rather it moves 
slowly with particles diffusing across the fluid-solid interface. The overall life-time of this 
crystal is several orders of magnitude larger than the external time-scale 7"^, imposed by 
the shear. The corresponding collision- angle distribution C{9) in Fig. 10(d) shows three 
distinct peaks at ^ = 0, ir/S and 27r/3. Note that the peak at ^ = n/S does not exist 
in the fluid phase (refer to Figs. 8a), and this provides evidence that the particles in the 
crystalline-phase are arranged in the triangular-packing structure. The large area of the 
crystaUine phase as well as its relatively large value of solids fraction (hence, small mean 
free time) further suggest that most of the coUisions occurred in the crystalhne-phase. 

We should remark here that, at such high densities, the inelastic collapse^^^"^^ would 
eventually terminate the evolution of the system. We used the TC-model"^*^ to avoid inelastic 
collapse, but could not altogether eliminate it within the crystal after some time. But the 
important point to note is that the crystalline-phase can be maintained for a long period of 
time (t > 100 strain units), and hence the reported results are not transient effects. 

Analogous plots for a nearly elastic system (e = 0.99) are shown in Fig. 11 (a-c) for the 
same density u = 0.8. The overall features are similar to that for e = 0.9, but the width of 
the crystalline zone is a little larger. By decreasing the dissipation-level to e = 0.7, we did 
not observe crystal formation, with other parameters being fixed; by increasing the system- 
size to = 4096, however, we observed layered crystal at e = 0.7. Thus, the formation of 
such layered-crystalline structure depends crucially on the system size and the dissipation 
level: the larger the system-size and/or the weaker the dissipation, the more susceptible the 
system is to crystallize. 

Note that even if we arc well below the limit of perfect square-packing (ugp = 7r/6 ~ 
0.785), the system could crystallize if the dissipation-levels are low; for example, we ob- 
served layered crystalline structures ai v = 0.75 with e = 0.99 and A^ = 4096. Decreasing 
the coefficient of restitution to e = 0.9, the flow-field remained homogeneous. Thus, our 
layered crystalline structures appear to be tied to a long-wave instability of the clastic 
hard-sphere fluids^^. Since our results are not driven by the inelastic dissipation, they are 
distinctly different from the layered shear-banding patterns, as predicted by the kinetic the- 
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ory models^^' ' , in gramilar CoTictte flows. Having said that we need to mention that such 
dissipative layering patterns (i.e. those which become stronger with increasing dissipation 
levels^^) were also found in simulations of a dilute sheared granular fluid"^^, but the simula- 
tions were allowed to evolve from an unsheared initial configuration with a special kind of 
boundary condition. In contrast, all our simulations started from a uniform shear condition. 
In the present contribution, we have mainly focussed on the non-dissipative layering in the 
dense limit, and the related issues of dissipative layering are relegated to a future study. 

Before moving further, we make a qualitative comparison with the earlier simulation 
work of Campbell & Brennen^^ on the bounded Couette fiow of inelastic, frictional hard 
disks. They also reported similar layered-microstructure but due to the small system-size 
{N = 40) and boundary-effects (they considered a shear-flow bounded by frictional walls), 
the distinct shear-band formation that we have reported is not evident in the snapshots of 
their simulations. Nevertheless, we believe that our results are akin to that reported by 
Campbell & Brennen. One of the referees has drawn our attention to the recent work of 
Campbcll"^^ who probed the dense limit of a three-dimensional Couette shear flow using soft- 
sphere simulations. He was able to maintain uniform shearing even at a density of u = 0.62, 
but beyond that he reported shear-band formation in that the flow-field degenerates into 
sheared and non-sheared zones in the gradient direction. This simply suggests that one 
can maintain uniform shearing in a three-dimensional geometry even beyond the analogous 
square-packing limit (0 = 7r/4 ^ 0.52) since the particles have now an additional degree 
of freedom, orthogonal to the shear plane, to rearrange themselves. Moreover, most of 
Campbell's simulations were done with 1000 particles (equivalently, 100 particles in 2D) at 
a restitution coefficient of 0.7. With these parameter vahics. wc did not find layered crystals 
at u — 0.8 in two-dimensions (see discussion in the next section). 

1. Shear-induced ordering and Reynold's dilatancy 

One of the referees has drawn our attention to the recent work of Lutsko^^ who studied 
the shear-induced ordering in a low-density {u 0.26) elastic hard-sphere fiuid. The earliest 
simulations of Erpcnbcck^^ on 3D clastic hard-sphere fluids showed that at high shear rates 
the system breaks down into orderded (solid) and disordered (fluid) phases in the direction 
of the mean vorticity (i.e. normal to the shear-plane). This induces a long-range two- 
dimensional ordering, called a string-phase, in the shear-plane. More importantly, such 
ordering occurs only for a range of shear-rates- the lower the density, the larger this shear- 
rate interval. However, the latter work of Evans and Morriss^^ showed that the disorder-order 
transition of Erpenbeck arises due to the profile-biased-thermostat^^ since the string-phase 
vanishes completely when a profile-unbiased-thermostat is used. Thus, the shear-induced- 
ordering in elastic hard-sphere fiuids appears to depend on the choice of the thermostat. 
Note further that all the above works probed moderately dense systems only, well below the 
corresponding square-packing limit. Now to compare these results with our observations in 
granular fluids, we flrst need to deflne an equivalent shear rate since the granular energy and 
the shear rate are dependent on each other for the latter system. Prom the energy balance 
equation, it is trivial to show that the granular energy has the following functional relation 
with the shear rate and the restitution coefficient: 

oc oc 7 , 

1 — 
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where 7* = 7/(1 — e^)^/^ is defined as the reduced shear rate. In the quasielastic hmit 
(e — 1) the reduced shear rate (and hence the granular energy) approaches infinity, which 
is equivalent to the high shear-rate limit of an elastic hard-sphere fluid. Since we observed 
layered crystalline structures only in the quasielastic limit, we may thus conclude that such 
structures would also persist in elastic hard-disk fluids at large shear rates. 

It is interesting to ask whether the layered crystalline-structures of Figs. 10 and 11 
are, in any way, related to the concept of the Reynold's dilatancy^'''^^ which is explained 
schematically in Fig. 12. The top two sketches depict the classical constant-load^^ shear-cell 
experiment (in two-dimensions) in which the material is subjected to a constant normal load. 
Fig. 12(a) corresponds to an ideal situation of maximum packing, with the particles being 
arranged in a triangular lattice (u = Utp = ti/2\/3). Clearly, in this situation the top and 
the bottom plates of the shear cell will simply slide over the material, without deforming 
it. However, if one of the plates is allowed to move in the vertical direction and thereby 
allowing the particles to rearrange themselves, e.g. as in Fig. 12(6), the material can be 
deformed even homogeneously (if ly < i^sp = 7r/6). This is the shear-coupled volume change, 
commonly known as the Reynold's dilatancy. Note, however, that our simulations mimic 
constant-volume experiments, since the volume of the computational box remains fixed (see 
Figs. 12c and 12d); but, of course, now the pressure (i.e. normal load) can vary in response 
to particle motions inside the shear cell. For this case also, the uniform shearing is possible 
if and only if the overall density remains below the square-packing limit (z/ < Usp). However, 
for densities above this value {usp < u < utp), the shearing can be started only if we allow the 
particles to rearrange themselves. This is possible if a part of the system becomes denser, 
allowing free volumes to the rest of the system which is nothing but Reynold's dilatancy 
too. Hence, we will end up with a crystalline-phase coexisting with a fiuid phase, a typical 
example of which is shown schematically in Fig. 12(d). Thus the phenomenon of Reynold's 
dilatancy, for densities Vgp < v < Vtp, would make the ordering transition, as depicted in 
Figs. 10 and 11, more prominent, and this effect would be much stronger in two-dimensions. 

IV. CONSEQUENCES FOR THE CONSTITUTIVE MODELING: 

RELAXATION MODELS 

Here we attempt to describe the normal-stress behaviour of a granular fiuid using the stan- 
dard relaxation-type models. Prior literature on the dense-gas kinetic theory, which forms 
the foundation of theoretical developments of granular fiuids in the rapid-shear regime, in- 
dicates that such a stress relaxation mechanism docs also exist in granular fiuids*^^"^^. The 
relaxation-type models arc routinely used to describe the non-Newtonian behaviour of vis- 
coelastic/viscoplastic materials, and hence might be apt for granular fiuids in the dense limit 
as well. The recent work of Zhang & Rauenzahn^^"^^ clearly shows that such viscoelastic 
stress relaxation mechanism exists in dense granular fiows. Following a rigorous statistical 
mechanical procedure, they derived an evolution equation for the collisional stress tensor 
which boils down to a frame-indifferent viscoelastic model, with the Jaumann derivative 
appearing directly without appealing to objectivity arguments. 

Let us consider the viscoelastic relaxation approximation suggested recently by Jin & 
Slemrod^^ to regularize the Burnett order equations of Sela & Goldhirsch^ for a low-density 
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granular fluid. Their proposed equation for the pressure deviator, in our notation, is 



n + n 



( 



DU 



Dt 



L^-n-n-L + |tr(n-L)iJ 
+ T2 (s • n + n • s - |ir(n -8)1)= n' 



eg 



(11) 



where 



s 



-2/iS - (AV-u) 1 + n2 + Hg, 
|(L + L-)-i(V.u)l, 



L 



A 




n 



0.3211 (^) . 




Here ri and T2 are relaxation times, d the dimensionality of the system, L is the velocity 
gradient, S the deviatoric part of the rate of strain tensor, /i the shear viscosity, ( the bulk 
viscosity and 1 the identity tensor; 112 said are higher order terms, explicitly written 
down in Jin & Slemrod^^. Note that both relaxation times are proportional to the ratio of 
the shear viscosity and the pressure, and hence proportional to the mean free time. In the 
limits of Ti, T2 — > and 112, ^3 — > 0, we recover the standard Newtonian model for the 
stress tensor. 

Neglecting the higher-order terms, an expression for the first normal stress difference can 
be obtained for the steady uniform shear flow: 



This quantity is always positive, as in our simulation results for dilute flows, if the two 
relaxation times are of the same order. 

It is important to note that the above evolution equation does not satisfy the principle of 
material frame indifference (MFI) which states that the constitutive laws should be invariant 
under rigid-rotation^'^. The scalar field 0, the vector field v and the tensor field II are called 
frame-indifferent or objective if the following relations hold for all t: 



where 0, v & 11 and 0', v' & 11' are defined in two different frames and respectively, 
and Q(t) is a proper orthogonal tensor. Here T' G S{J-'), with S{J-') denoting the set of all 
frames obtainable from a given frame by observer transformations. That the stress-tensor 
in a granular gas is not a frame- independent quantity (as in the hard-sphere gas^"^'^^) is well- 
known. Since the kinetic component of the first normal stress difference remains positive 




> 0. 



(12) 




Q(t)v(x,t), 
Q(i)n(x,i)Q(i)^, 



(13) 
(14) 
(15) 
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at all densities, the kinetic stress tensor can be modelled using a non-objective equation 
as discussed above. For the collisional stress tensor, one can postulate a similar evolution 
equation as in Eqn. (11), but ti must be multiplied by a factor which must change sign 
at the critical density. However, this equation would remain frame-dependent even in the 
dense limit. Thus, if one has to recover the Boltzmann-limit of relaxation- type equations, a 
frame-indifferent approach does not appear to work. 

A question now arises as to the possibility of modelling normal stress differences using 
the standard frame-indifferent relaxation type models. The simplest way is to use either the 
lower-convected or the upper- convected equations for the pressure deviator: 

n + a{iy, e)T + {L^ • H + H • L - |tr(n • L)l}^ = -2/xS - (AV-u) 1, (16) 

n + a(i/, e)T - {L • n + n • L'^ - |tr(n • L)l}^ = - (AV-u) 1, (17) 

respectively. Here q;(i/, e) is an empirical constant, dependent on both the density and 
restitution coefficient. For both cases, the first normal stress difference is 



1 -I- r^a^ \p 



Clearly, if q;(z/, e) is obtained from simulation, its sign-reversal would also correspond to 
the sign-reversal of jVi. Comparing the Jin-Slemrod equation with its corresponding frame- 
indifferent analog (lower-convected model), we conclude that the loss of frame-indifference 
shows up as a sign- change of the first normal stress difference. It would be interesting to 
investigate whether one could relax the Sela-Goldhirsch equations^ using a frame-indifferent 
approach without violating the entropy inequality^^ . 

Similarly, one could postulate evolution equations using other objective derivatives. In 
this regard, the co-rotational Jeffrey's model seems to be the ideal choice: 

n + n(z., e)^ = -2/. + - (V-u) 1 + r,{u, e)— j (19) 

with V/Vt being the Jaumann derivative^. The corresponding first normal stress difference 
is 

M^-%^(^:). (20) 

which is positive/negative depending on whether ti is less/greater than T2. Thus, the frame- 
indifferent relaxation models are able to predict positive and negative first normal stress 
differences. (For the steady homogeneous shear flow, one can also model positive/negative 
normal stress differences by postulating a general orthonormal basis, generated by the nilpo- 
tent basis tensors, which satisfies the objectivity requirement; for related issues, the reader 
is referred to Goddard^.) 



V. SUMMARY AND CONCLUSION 

We have probed the non-Newtonian behaviour and the incipient crystalline-phase in a 
sheared, monodisperse, two-dimensional granular fluid. The standard event-driven tech- 
nique is used to simulate a box of hard-disks under homogeneous shear deformation. The 
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information about the stress tensor is obtained by decomposing it in the standard way: 
P = pi + H, where p is the pressure and 11 the pressure deviator. The non-Newtonian 
behaviour is quantified in terms of the first normal stress difference: Mi = {Hxx — ^yy)/P- 

The granular fluid is non-Newtonian with a measurable flrst normal stress difference 
(A/i) which is positive (if the stress is defined in the compressive sense) in the dilute limit. 
Interestingly, however, A/i changes from positive to negative at a critical density in the dense 
regime. By decomposing J\fi into the kinetic and coUisional contributions, A/i = AAf + Afi, 
we found that while J\f-^ is always positive and decays to zero in the dense hmit, has a 
non-monotonic variation with density. In particular, jVf increases from zero in the dilute 
limit as v increases, reaches a maximum at some value of v and then decreases, eventually 
becoming negative in the dense limit. The density at which A/'f = (z/ = 77) depends 
crucially on the level of micro-scale dissipation; in particular, V increases as the coefficient of 
restitution decreases. We have constructed a phase-diagram in the (z/, e)-plane by identifying 
the regions where Mi is positive/negative. 

We have shown that the origin of the first normal stress difference, in the dense limit, 
is tied to shear-induced coUisional anisotropics. The underlying mechanism is distinctly 
different from that is known for a dilute granular gas^'^^ where the anisotropy in the second 
moment of the fluctuation velcoity, which is a Burnett-order effect, gives rise to normal 
stress differences. 

At the micro-level, the particles undergo reorganization as the dense-limit is ap- 
proached. The signatures of microstructural-reorganization have been captured by probing 
the collision- angle distribution, C(^), which is anisotropic due to the presence of the mean 
shear. In particular, we have found that the topology of the anisotropic-structure of C {6) 
changes, with collisions occurring at certain preferred angles on the upstream-faccs of the 
colliding pairs. The sign-reversal of Mi is correlated with a preferred value of the average 
collision angle, Oav = 7r/4 ± 7r/2, averaged over the upstream-ieices of the colliding particles. 

The time evolution of the sheared granular fluid leads to crystallization in the dense 
limit, signalling the coexistence of fluid and solid. The particles are arranged in a triangular- 
packing inside the crystal, and it moves as a layer in the strcamwisc direction. The formation 
of such layered-crystalline structure depends crucially on the system size and the dissipation 
level: the larger the system-size and the weaker the dissipation, the more susceptible the 
system is to crystallize. This appears to be related to a long-wave instability^^ of the elastic 
hard-sphere fluids. 

The present work clearly shows that the available kinetic-theory-based rheological models 
for granular fluids are not adequate to predict the behaviour of the flrst normal stress 
difference in the dense limit. Certain microstructural-features, like the preferred distribution 
of collisions which eventually leads to crystal-formation, should be incorporated into the 
theory. At such high densities, many-body effects (both positional and velocity correlations) 
are important^^'"^^ and the appropriate kinetic description is the BBGKY-hierarchy^^. To 
incorporate the observed normal stress behaviour into the framework of plausible constitutive 
models, we showed that the standard frame-indifferent relaxation type models can be used 
to model both positive and negative first normal stress differences. In this regard, the two- 
parameter Jeffrey's modef appears to be the ideal choice; however, we are unable to recover 
the corresponding Boltzmann limit which is known to be non-objective. On the whole, we 
believe that a lot remains to be done for a better understanding of the dense-phase rheology 
of granular fluids even in the hard-sphere limit. 
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APPENDIX A: SHEAR-INDUCED ANISOTROPY AND THE COLLISION 

ANGLE DISTRIBUTION 

Here we are interested in the shear-induced anisotropy of the collision angle distribution 
of an inelastic hard-disk fluid. Following Savage & Jeffrey^^, an expression for the collision 
angle distribution C{9) is derived, which is compared with the simulation data in nB.2. 

To calculate the probabhhty of collisions at a specific angle 9, we focus on Fig. 13 with 
two particles colliding at r. Note that 6 is measured anticlockwise from the positive x-axis. 
For collisions to occur in a time 6t, the center of particle 1 must lie inside the volume 
(T5k(q • k)5t, where q = Ci — C2 is the relative velocity of the colliding pair. Thus the 
expected total number of collisions (per unit time and unit area) with the line of centres k 
lying between k — 5k/2 and k -|- 5k/2 is given by 

J af^\ci, ri, C2, r2)(q • k)(k • n)dkdcidc2, (Al) 

where f^'^\-) is the two-particle distribution function which is defined so that 
/*^^^(ci, ri, C2, r2) dcidc2dridr2 is the number of pairs of particles such that the particle 
i is located in an area element dri about ri with its velocity in the interval dci about Ci 
while particle j is located in an area element dr2 about r2 with its velocity in the interval 
dc2 about C2. To progress further, we have to invoke the assumption of molecular chaos and 
hence the expected number of collisions is 

J ag{ri, r2)f^'\ci, r,; u(ri))/«(c2, rs; u(r2))(q • k)(k • n)dkdcidc2, (A2) 

where g{ri, r2) is the pair-distribution function. For the steady uniform shear flow, g{ri, r2) 
is calculated from the relation^^'^^'^^: 

g{r„r2) = ^^ [ /«(ci,ri;u(ri))/«(c2,r2;u(r2))cicicic2, (A3) 

^ ^qk>0 

where gdi^) is the contact value of the pair-distribution function and q • k > implies that 
the integration be carried out for impending collisions. 

As a first approximation, the single particle velocity distribution function 
/*^^)(ci, ri; u(ri)) is assumed to have the Maxwellian-form: 

/(..(c..r.uW)^(^)exp 



m(ci - u(ri))' 
2kRT 



(A4) 
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where T is the granular temperature (fluctuating kinetic energy) and fc^ the usual Boltz- 
mann constant. Now transforming the particle velocities (ci, C2) to their center-of-mass 
and relative velocities, equation ()A3|) can be integrated to arrive at the following explicit 
expression for the pair-distribution function^^: 



fl'(ri,r2) = fi'c(z^)erfc 



2k-u(r2) 

{2kBTy/\ 



(A5) 



where erfc(-) is the complementary error function. Using the above expression for the pair- 
distribution function and transforming in terms of polar coordinates {r,6), the integral for 
the normalized collision angle distribution yields'^^ 



C{e) = AiT) 



exp 



sin^6'cos^^\ 5f(^) sin^cos^ 



2T / ^/T 
where g{6) is the angular pair-distribution function given by 



9(0), 



(A6) 



9(0) 



9[ri,r2) 
9c{j^) 



erfc 



sin 6 cos 6 



2T 



(A7) 



and A{T) is a normalization constant. 

For the uniform shear flow, an expression for the granular temperature, can be obtained 
from the energy balance equation, by equating the energy production due to shear-work 
with the energy loss due to colhsions: 



T 



(A8) 



where /i = f^{i')\/T is the shear viscosity and T> = (pp/cr)/x)(i/, e)T^/^ the collisional 
dissipation rate, with 



4 



i^9c{i') 



+ 2 + iyg,{u) 1 



TT 



Substituting this expression for T, the normalized collision distribution function becomes 

fv{^, e) sin^ 9 cos^ 9^ 



C{9; z/, e) = A{T) 



exp 



g{9) sin 6* cos 



9(0) (A9) 



and the angular pair-distribution function 



g{9; u, e) = erfc 



sin 6' cos 4 



(AlO) 
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It is clear that the angular dependence of C{9] i/, e) is modified by both the inelastic dissi- 
pation and the density. 
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FIG. 1: (a) A snapshot of the sheared granular system at steady-state. The arrows indicate the 
displacement of the image boxes. (6) Variations of the granular energy T and the calculated shear 
rate ^cai with time. For an explanation of the system and particle properties, see the text. The 
parameters for both subplots are v = 0.5, e = 0.7 and = 1024. 
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FIG. 2: (a) Variation of the first normal stress difference A/i with the solid volume fraction. The 
arrows on the left ordinate indicate corresponding analytical values for a two-dimensional granular 
gas. (b) Variations of A/i, Ni and AAf with at e = 0.7. In both subplots, the symbols represent 
the simulation data and the lines are drawn to guide the eye. 
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FIG. 3: The phase diagram, dehneating regions of positive and negative first normal stress differ- 
ence, in the (i/, e)-plane. The filled circles and triangles represent zeros of M\ and A/i, respectively. 
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FIG. 5: Radial distribution function g{r) plotted against the normalized distance r/a from sheared 
simulations with dissipation e = 0.7 (solid lines) and from homogeneous, non-sheared situations 
with 7 = and e = 1 (dotted lines): (a) u = 0.6, (6) u = 0.7; (c) u = 0.725; (d) u = 0.75. The 
arrows indicate the peak values of g{r) at contact for the homogeneous system. 
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FIG. 6: Variations of the pressure and viscosity functions with density at e = 0.7: fp = p/pT and 
ffj. = n/pfrVr. The hues are drawn to guide the eye. 
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FIG. 7: Distribution of collision angles C{9) for different coefficient of restitutions at = 0.6. The 
symbols represent simulation data and the lines theoretical predictions. 
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FIG. 9: Variations of the average collision angle dav with density for different restitution coeffi- 
cients. The arrows indicate densities where Mi ~ 0. The lines joining the data points are to guide 
the eye. 
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FIG. 10: Evidence of crystallization in a sheared dense granular fluid at = 0.8 and e = 0.9. (a) 
Particle distribution and (b) streamwise velocity at t = 390; (c) evolution of streamwise velocity 
at early times; (d) collision angle distribution. 
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FIG. 11: Effect of the coefficient of restitution on crystallization: p = 0.8 and e = 0.99. (a) 
Particle distribution and (b) streamwise velocity at t = 150; (c) collision angle distribution. 
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FIG. 12: A schematic diagram to explain Reynold's dilatancy in the Couette shear flow. 
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FIG. 13: A schematic of the possible collision angles 9 in uniform shear flow; 9 is measured 
anticlockwise from the positive x-axis. Note that the collisions are more likely to occur in the 
second {tt /2 < 9 < tt) and fourth (37r/2 < 6 < 27r) quadrants of the colliding disks. 
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